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We derive a coupled set of equations that describe the non-equilibrium evolution of cumulants of 
critical fluctuations for space-time trajectories on the cross-over side of the QCD phase diagram. 

In particular, novel expressions are obtained for the non-equilibrium evolution of non-Gaussian 
Skewness and Kurtosis cumulants. Utilizing a simple model of the space-time evolution of a heavy- 
ion collision, we demonstrate that, depending on the relaxation rate of critical fluctuations, Skewness 
and Kurtosis can differ significantly in magnitude as well as in sign from equilibrium expectations. 
Memory effects are important and shown to persist even for trajectories that skirt the edge of the 
critical regime. We use phenomenologically motivated parameterizations of freeze-out curves, and 
of the beam energy dependence of the net baryon chemical potential, to explore the implications of 
our model study for the critical point search in heavy-ion collisions. 


I. INTRODUCTION 

The structure of the QCD phase diagram has attracted 
a large number of theoretical and experimental stud¬ 
ies w- Of high interest is the possible existence of 
a conjectured critical point 01 a in the phase diagram. 
This critical point is the end point of a first-order phase 
transition line that separates, in the chiral limit, a chi- 
rally symmetric quark-gluon plasma (QGP) phase from 
a hadron matter phase of QCD. An entire experimental 
program, the Beam Energy Scan (BES) at the Relativis¬ 
tic Heavy-Ion Collider (RHIC), is dedicated to searches 
for the QCD critical point mm- 

A universal feature of a system near the critical point 
is the growth and divergence of the fluctuations of the 
order parameter. These fluctuations can be quantified 
by the variance of the critical field (the Gaussian cumu- 
lant), as well as higher non-Gaussian cumulants such as 
skewness and kurtosis, which correspond, respectively, to 
the third and fourth cumulants. They scale with the cor¬ 
relation length £, which is large near the critical regime 
and divergent at the critical point. 

As pointed out in Ref. [8], the non-Gaussian cumulants 
are much more sensitive to £ than the variance. For ex¬ 
ample, while the variance grows as £ 2 , the kurtosis grows 
far more rapidly as £ 7 . Further, even qualitative features 
of the non-Gaussian cumulants, such as a change in sign 
and the associated non-monotonicity, can signal the pres¬ 
ence of criticality in the QCD phase diagram HE]. As 
observed in Ref. HQ , the sign of kurtosis in equilibrium is 
negative when the critical point is approached from the 
crossover side and positive when approached from the 
first order transition side. 

These enhanced near critical fluctuations are accessible 


through measurements of event-by-event fluctuations of 
various conserved quantities in heavy-ion collisions m 
US] as well as fluctuations of particle multiplicities m 
eng. As the QCD phase diagram is scanned by varying the 
beam energy, non-Gaussian event-by-event fluctuations 
of multiplicities are generically expected to show non¬ 
monotonic behavior in the proximity of a critical point. 
Uncovering such behavior, and cleanly identifying this 
behavior as a signature of criticality, is a major focus of 
the RHIC BES program p2Hl9j. 

The above expectations are entirely based on the as¬ 
sumption that the soft modes responsible for critical fluc¬ 
tuations are in equilibrium with the medium. Indeed 
if this were the case, fluctuations of conserved quanti¬ 
ties measured in heavy-ion experiments could be directly 
compared to equilibrium lattice QCD calculations |20f- 
[23] . However the expanding medium created in heavy- 
ion collisions only spends a limited amount of time in the 
QCD critical regime and it is unlikely the critical modes 
remain in equilibrium in this duration. Therefore the rel¬ 
evant cumulants can in principle differ considerably from 
their equilibrium values mis]. For instance, the time it 
takes the correlation length to reach its equilibrium value 
scales as r e ff ^ £ 2 , which defines the dynamic scaling ex¬ 
ponent £ [23 BE] • 

Since the non-Gaussian cumulants grow with higher 
powers of the correlation length, their relaxation times to 
their equilibrium values may be significantly larger. Even 
the sign of a non-Gaussian cumulant may differ from its 
equilibrium value due to this critical slowing down of re¬ 
laxation rates. These simple considerations suggest that 
non-equilibrium memory effects may play an important 
role in interpreting the results of BES measurements, and 
conversely, in applying lattice results to predict detailed 
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features of the data. 

A number of papers have previously investigated crit¬ 
ical dynamics off-equilibrium. The quasi-stationary dy¬ 
namics of Gaussian fluctuations near the critical point 
was studied in Ref. [24] by using a model equation to 
compute the relaxation of the correlation length. Gaus¬ 
sian fluctuations in presence of the critical point were 
also investigated within a hydrodynamic approach m- 
Further, the hydrodynamic evolution of a quark-gluon 
fluid that is coupled to the classical real-time evolution of 
long wavelength modes of chiral fields has been explored 
in a number of effective models that contain a critical 
point [28H3Cf . These models have also been employed to 
investigate the out-of-equilibrium evolution of the rele¬ 
vant correlation length. Another class of models employs 
a coupled Boltzmann-Langevin-Vlasov kinetic approach, 
with conserved charges to model the evolution of critical 
fluctuations through chemical and kinetic freeze out m • 
These can be contrasted [32 with results from UrQMD, 
a hadronic event generator that does not contain criti¬ 
cal fluctuations, to isolate genuine effects due to critical 
fluctuations. All the studies outlined focus on the evo¬ 
lution of Gaussian cumulants. (See Refs. |33j [34] for the 
evolution of higher order cumulants in a system with¬ 
out critical phenomena.) Analyses incorporating critical 
fluctuations of non-Gaussian cumulants off-equilibrium 
are sorely lacking. 

This paper is a first attempt to address this gap in 
our knowledge of non-Gaussian fluctuations and to un¬ 
derstand its consequences. In particular, we compute 
the real-time evolution of non-Gaussian cumulants of the 
critical field to estimate memory effects in the vicinity of 
the QCD critical point. The expressions for the cumu¬ 
lants are derived from the Fokker-Planck equation deter¬ 
mining the evolution of the non-equilibrium probability 
distribution function of the critical field a. One obtains 
a set of coupled first order differential equations which 
describe the evolution of the mean, variance, skewness 
and kurtosis, on the cross-over side of the QCD critical 
scaling regime, for trajectories that are close to but not 
at the critical point. This coupled set of equations takes 
into account universal equilibrium properties and non¬ 
equilibrium dynamics in the QCD critical regime and is 
the appropriate theoretical framework to study the evo¬ 
lution of non-equilibrium cumulants in that regime. 

The QCD critical point lies in the static universality 
class of the 3-dimensional Ising model [5] for which the 
equilibrium cumulants are well known [35]. The latter 
can therefore be used to fix the equilibrium values of the 
cumulants of the QCD critical field. With these fixed, 
the non-equilibrium values of the cumulants can be de¬ 
termined from our evolution equations. A key quantity 
determining the non-equilibrium evolution of cumulants 
is the effective relaxation time r e ff of critical fluctuations. 
This quantity is unknown and can only be estimated to 
be a strong interaction time scale. For reasonable varia¬ 
tions in its value, we will show that strong memory effects 
are seen for the non-Gaussian cumulants. 


While our expressions in terms of Ising variables are 
universal, the relation of these to the temperature T and 
the chemical potential /i in QCD is non-universal, and 
is a source of significant systematic uncertainty. Never¬ 
theless, with physically motivated assumptions, one can 
obtain expressions for the cumulants as a function of tem¬ 
perature and chemical potential. There is a significant 
amount of phenomenological information from thermal 
model ratios of particle multiplicities that allows one to 
extract the chemical freeze-out curve in temperature and 
chemical potential in the QCD phase diagram. Further, 
the particle ratios, with model assumptions, can be used 
to relate fi to the center of mass energy yfs of the heavy 
ion collision. With these phenomenological inputs, our 
results can be used to provide qualitative estimates of 
the importance of memory effects for the BES. 

We will demonstrate that, depending on the trajectory 
followed, the sign of Skewness and Kurtosis can flip rel¬ 
ative to their equilibrium values. These results suggest 
the need for caution in interpreting, on the basis of equi¬ 
librium expectations, the results of experiments. We also 
note that if the system were in equilibrium throughout 
its evolution in the T-/i plane, it would have to pass very 
close to critical point to be sensitive to critical fluctu¬ 
ations. On the other hand, if memory effects are im¬ 
portant, even a trajectory some distance away from the 
critical point may be sensitive to it. 

The paper is organized as follows. In Sec. [II] we review 
the equilibrium properties of the a field in the critical 
regime and derive the non-equilibrium evolution equa¬ 
tions for cumulants. Our results up to this point are 
completely general and apply equally to all systems that 
lie in the 3-D Ising universality class. In Sec. Ill A, we 


map our results to the QCD critical regime. Because 
the QCD critical point lies in the universality class of 
the 3-D Ising model, we will be able to specify the equi¬ 
librium properties that are necessary inputs for apply¬ 
ing our equations to study the evolution of cumulants in 
the QCD critical regime. We construct a simple model 
of the space-time evolution of a heavy-ion collision in 


Sec. |IIIB| In Sec. |III C[ we solve the evolution equations 
and present results for a representative trajectory pass¬ 
ing through the critical regime. We further extend our 
formalism to model the beam energy scan in Sec. |IV[ We 
summarize our results and discuss their implications for 
the QCD critical point search in Sec. [V] 


II. EVOLUTION EQUATIONS FOR 
CUMULANTS 


A. Critical fluctuations in equilibrium 


We begin by considering the critical field cr(x) and its 
zero momentum mode 




t J d 3 xa(x), 


( 2 . 1 ) 
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where V is the volume. The fluctuations of this zero¬ 
mode cr-field are described by the probability distribution 
P(cr; r) which in general depends on the (proper) time 
r. Introducing the average of any cr-dependent functions 
with respect to P(cr;r), 


f-oo da ' (•••) P ( <J ’ T ) 

JZo daP ( a ; r ) 


( 2 . 2 ) 


of <j(x). Consequently, the equilibrium distribution be¬ 
comes 

P 0 (a) ~ exp (-VMa)) , V A = ^ , (2.6) 

where Qo(,ff) is a function of the zero mode of the cr-field, 
^o(c) = ^m^(cr-cro) 2 +y(o'-(To) 3 +h((j-(Jo) 4 - (2.7) 


one can define the expectation value of the cr-field M(the 
“magnetization”) and the cumulants of the cr-field as 

M(r) m (a ), Scr = a — M(r) 
k 2 (t) = ({5a) 2 ) , k 3 (t) = ({5a ) 3 ), 

k a {t) = ((5a) 4 ) - 3 [k 2 (t)} 2 . (2.3) 

Before discussing the evolution of non-equilibrium crit¬ 
ical fluctuations, it is useful to first review their proper¬ 
ties in equilibrium. For the static universality class of 
the 3-D Ising model, the dependence of the equilibrium 
cumulants on the Ising variables (the reduced tem¬ 
perature r = (t — t c )/t c , where t is the Ising temperature, 
t c the Ising critical temperature and h the rescaled mag¬ 
netic field) is universal. The explicit expressions for these 
are given in Appendix |A| 

Alternatively, the critical fluctuations in equilib¬ 
rium can be described by the distribution P 0 (cr) ^ 
exp(—£o(a)/T), where T is the temperature. Here we 
expand the effective action functional £o(a) around the 
minimum cr 0 in gradients and powers of cr(x) as 

£o(cr) = J d 3 xl i [Vcr(a :)] 2 + t to 2 [a(x) - a 0 ] 2 + 

y [a(x) - <T 0 ] 3 + y K*) - Oo ] 4 j > 

(2.4) 

following Refs. mm- Near a critical point, the mass 
(m a ) of the cr field, as well as the other parameters in 
Eq. ( |2.4| ), scale with the equilibrium value of the correla¬ 
tion length £ eq as 


An important quantity characterizing the probability 
distribution Eq. (2.6) is the ratio between the correlation 


length and the size of the system, 

fW 


(2.8) 


Throughout this work, we will work in the scaling regime 
(near, but not at the critical point) where the correlation 
length £ is larger than any microscopic scale P m i cr but 
smaller than the size of the system, P m i cr £ < C L . 
Hence e < 1. In this regime, the kinetic term in Eq. (2.4) 
is proportional to cr 2 /P 2 and is small compared to the 
mass term cr 2 /£ 2 n . This justifies our dropping the kinetic 
term in Eq. (2.4). Further, truncating the expansion at 


the A 4 term in Eq. (2.4) is justified as well because higher 


order terms are suppressed for small e. 

Treating e as an expansion parameter enables us to 
relate the parameters ao,m a , A 3 , A 4 to the equilibrium 
“magnetization” M eq and the cumulants ft® q ,n = 2,3,4. 
To leading power in e, we obtain, 


M eq = cr 0 , 


/4 q = 


eq 


£ 2 

Seq 

V 4 

6 


eq 

£* ^ - 


2 A, 


= 


6 

eq ' 


y3 [2(^3£eq) ^ 4 ] £ e q 5 (2-9) 


where we have employed Eq. (2.2), ( |2.3| ) and (2.7). This 
result is in agreement with Ref. [ 8 f Using Eq. (2.5), 
we also recover the scaling behavior M eq 


_1 / 2 „eq 


£ 2 

Se^ 


6 


. 9/2 


Seq 


B. Evolution equation for cumulants 


m cr 1 = £eq , <?0 = CT 0 T(T£ eq ) 1/2 , 

A 3 = A 3 T(T £ eq )- 3 / 2 , A 4 = MT&d)- 1 ■ (2.5) 

Expressed thus, the dimensionless parameters do,A 3 ,A 4 
do not depend on £ eq and therefore remain finite when 
approaching the critical point. Because in this work we 
are interested in the fluctuations of the zero momentum 
mode of the cr field, we will neglect the spatial dependence 


1 Following Refs. min], we shall neglect the anomalous scaling 

dimension 77, which is only of order few percent, in this work. 


Now that we have established the behavior of the equi¬ 
librium cumulants in our approach, we will turn to their 
non-equilibrium evolution with proper time r. The non¬ 
equilibrium transport properties of fluids near a critical 
point depend on their dynamical universality class; these 
were formalized in the classification scheme of Ref. [S5j . 
In the QCD case, which according to (26} f36il38l should 
lie in the model H universality class of [25] . the relevant 
fields are those of the chiral order parameter (the chiral 
condensate) and the baryon density. 

It was shown in [26] that the relevant hydro dynam¬ 
ical equations for the space-time evolution of these two 
fields can be expressed as a coupled set of Langevin equa¬ 
tions. In the late time, long wavelength hydrodynamic 
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asymptotic limit, the corresponding eigenvalue problem 
can be solved, revealing only one long wavelength diffu¬ 
sive mode, with a diffusion constant that goes to zero at 
the critical point. As argued in Ref. (26] , both the chiral 
condensate and the baryon density fluctuate on these hy¬ 
drodynamic time scales, with the chiral condensate seen 
as “tracing” the profile of the baryon density as it relaxes 
to its equilibrium value. 

The diffusive properties of this critical mode are cap¬ 
tured by the Fokker-Planck equation for the relaxation 
to equilibrium of the distribution P(cr, r), expressed as 


f^o(cr) as a function of a around a = M: 

n' 0 (a) = Q' 0 (M) + n%(M)5<r 
Ua' 


W{M)_ ga3 


2 ! 


W(M). 3 


3! 


(2.13) 


As Qo(a) is a polynomi al of cr, the ab ove ex pansion is 
exac t. S ubstituting Eq. (2.13) into Eq. (2.12) and using 
Eq. (|2.3|), we obtain, 


<9 r M(r) = - 


1 


d T P(a ; t) = 


1 


(m^Teff) 


d a [9 <T n 0 (o') + V 4 1 da] P(cr, t) 






( 2 . 10 ) 


3! 


k 3 (t) 


(2.14) 


where r e ff is the effective relaxation rate. As we shall 
discuss further shortly, r e ff scales with a universal power 
of the correlation length £ that is characteristic of model 

H. 

In a static medium where T an d V are independent of 
time, it is e asy to check that Eq. ( |2.6| ) is the static solu¬ 
tion to Eq. (2.10). For an expanding medium, the solu¬ 


tion is more involved. Now Qo(a) also depends on time 
since the thermodynamical variables such as the temper¬ 
ature and the chemical potential change with time. Nev¬ 
ertheless, as long as the expansion rate does not exceed 
the characteristic time scales in the system, Eq. (2.10) 


will describe the evolution of the distribution towards a 
quasi-stationary fixed point solution. 

We shall now derive the evolution equations for the cu- 
mulants. We first note that for any function g(a) which 
does not grow exponentially in the large a limit, the evo¬ 
lution of the expectation value ( g(a )) is given by 




/ oo 

dag(a)d T P(a ; - 

-OO 

(/: 


■) 


dag{a) [SI' 0 (<j)P(<j; r)] 

) 

/ OO 

dag(a)P"(a-,T) 

-OO 

<5»>f 


(m^Teff) . 


(g'(a)fl' 0 (cr)) - 


V 4 


, ( 2 - 11 ) 


where we have used Eq. (2.10) and the definition 


Eq. ( 2 . 2 ). Here and hereafter, we will use the prime 
symbol to denote the derivative with respect to cr. To 
obtain the last line in the equality above, we performed 
an integration by parts. 

We will firs t stud y the evolution of M by taking g(a) = 
a. From Eq. (2.11), we immediately obtain 


<9 r M(r) = - 


1 


( m a T es) 




( 2 . 12 ) 


whe re we have used Eq. ( |2.3| ). To express the R.H.S of 
Eq. (2.12) in terms of cumulants, we first Taylor expand 


The evolution equation Eq. (2.14) up to this point is ex¬ 


act but not closed as the R.H.S of Eq. (2.14) depends 


on the non-equilibrium values of , ft 3 . We shall now 
demonstrate that the the evolution of ^ 2,^3 decouples 
from the evolution of M(r) for small values of e. 

In order to arrive at this result, we first introduce the 
dimensionless functions 


F n (M) = V 4 [e 2 ~ n b n d^ 0 (<t)] 


<j=M 


n = 1,2,3,4, 
(2.15) 


whic h, by construction, are finite in the small e limit. In 
Eq. (2.15), the dimensionful quantity b is the square root 


of the variance in Eq. (2.9) 



(2.16) 


The reader should keep in mind that that, for an ex¬ 
panding medium, b will depend on r due to the change 
of both V and £ e q with proper time. Explicitly working 
out Eq. ( |2.15| ), we obtain 

F 1 (M) = 6M [l + A s (SM) + A 4 (£M) 2 

F 2 (M) = 1 + 2 A s (SM) + 3A 4 (£M ) 2 , 

F 3 (M) = 2 |~A 3 + 3A 4 (£M)1 , F 4 = 6 A 4 . (2.17) 

On the R.H.S, we have defined, for convenience, the di¬ 
mensionless q uant ity SM = e (M — a 0 ) /b. 

Using Eq. (2.3), the evolution equation for M(r) can 
be expressed as 

drM( t) =-t- s 1 (('){f 1 (M) 


k 18 > 


We shall now consider Eq. (2.18) in the small e limit. 
In powers of e, we will count 

eM k 2 k 3 «4 nm / 91Q x 
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This power counting, as is clear from Eq. (2.9), holds for 
the equilibrium cnmn lants. In the following, we will as¬ 


sume that Eq. (2.19) also holds for the non-equilibrium 
cumulants. To the extent that the Fokker-Planck “mas¬ 
ter” equation is valid, this assumption is reasonable. We 
will confirm later that this power counting is consistent 
with our evolution equations at all times. Consequently, 
the and terms on the R.H.S of Eq. (2.18) are sup- 
e 2 

closed form expression for M to be 


the Gaussian distribution, 

1 


tto(cr) = -m a (a - <r 0 ) • 


( 2 . 21 ) 


For this case, the equilibrium values of K n are simply, 


M eq = (TO , 


= b 2 


eq 


/C = o. 


( 2 . 22 ) 


pressed by e 2 and e 4 respectively and one obtains the 


Further, Fi = 5M,F 2 = 1, F 3 = F 4 = 0, and Eqs. (2.20) 
reduce to 


d T M(r) = -t~s Fi(M) [1 + 0(e 2 )] .(2.20a) d r K n = -nrj [k„(t) - < q ] , n = 1,... (2.23) 


We wish to emphasize again that while 6, e is independent 
of time for a static medium, for the expanding medium 
we will study in this paper e, b will be time dependent due 
to the change of volume V(r), temperature T(r) and £ eq . 

The above derivation can be straightforwardly ex¬ 
tended to obtain evolution equations for higher cumu¬ 
lants. The details are given in Appendix. [B] In general, 
the evolution of K n is coupled to both lower cumulants 
hv n —i , hv n —2 and higher cumulants /c n+ i,/c n+ 2 . However 
the coupling to higher cumulants is suppressed by pow¬ 
ers o f 6. Keeping contributions at leading order in e from 
Eq. (B8) in Appendix. [B| we obtain the expressions, 


cWr) = —2r-ff 1 (b 2 ) [(^f) F 2 (M) - i] [1 + 0(e 2 )] , 

<Wr) = -3^ (, sb 3 ) [(^|)F 2 (M)+ (^) 2 F 3 (M) 
x [1 + d( e 2 )] , 


a T K 4 (r) = —4 r 



+ 3 (?) (SO Fs(M) + 

x [l + 0(e 2 )] . 


(S)‘ 


F 4 


(2.20b) 


As the R.H.S of Eqs. Q2.20 ) only depends on M, /t n , n = 
2, 3,4, this system of equations is closed. Eqs. (2.20) are a 


key result of this paper, and to best of our knowledge, are 
new in the literature. Employing the power counting in 
Eq. (2.19), one observes that the R.H.S of the evolution 
equations satisfies this power counting. Therefore if ini¬ 


tially the power counting Eq. (2.19) is satisfied, it is pre¬ 
served for all subsequent times. Since in our derivation 


we do not assume that the medium is static, Eqs. (2.20) 


are well suited for studying dynamical systems such as 
those created in heavy-ion collisions. 


C. Two limiting cases 

It is in structive to examine the evolution equations 
Eq. (2.20) in limiting cases. We first consider the limit 


where the equilibrium probability distribution Po(cr) is 


where note that ki is shorthand for the magnetization 
M. Eq. (2.23) expresses the fact that if the equilibrium 


probability distribution of the cr-fleld is a Gaussian, the 
evolution of cumulants are decoupled. For ft n , the damp¬ 
ing rate is proportional to n; hence the higher cumulants 
approach their equilibrium values earlier than lower cu¬ 
mulants. 

With certain assumptions, the results in this Gaus¬ 
sian limit can be shown to be identical to those obtained 
previously in Ref. m- The latter follow from the rate 
equation conjectured to describe the evolution of the non¬ 
equilibrium correlation length, 


9t [£ V)] = -tJ [£ X (r) -£eqV)] • 


(2.24) 


To facilitate a comparison of this equation to Eq. ( |2.23| ) , 
one identifies the non-equilibrium correlation length to 
be 


£(t) = vV 4 k 2 (r) ■ 


Eq. (2.23), for n = 2, can then be expressed as 

'e(ry 


d T 


V 4 


= -2 T 


-1 

eff 


v 4 


f 2 

Seq 


(2.25) 


(2.26) 


If we require the medium to be static and require further 
that deviations of £(r) from t he eq uilibrium value ^ eq are 
small, we then find that Eq. ( 2.26| ) reduces to Eq. (2.24) 
after a rescaling of r e ff by a factor of 2. 

Another interesting limit is the near equilibrium limit 
where Sk u — K n — becomes small. In this case, the 
evolution equations can be linearized to read, 


Rff 1 




d r M(r) = 

d T ^{r) = -3 r e y (e& 3 ) ( ) + 4A 3 ( ^ 


d T MT)] = -4T- ff V& 4 )< 


d T K 2 (r) = -2r eff 1 (5K2 


- 6 (2A| - 3A 4 ) |. (2.27) 


From the power counting in Eq. ( 2.19), all the terms in 
the [...] of the R.H.S of Eq. (2.27) are in the same order 
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in power of e. One can make the following observations 
about the evolution of cumulants in this particular limit¬ 
ing case. Firstly, unlike the Gaussian limit, the evolution 
of higher cumulants are coupled to the evolution of k 2 - 
It follows further, in contrast to the Gaussian limit, that 
the higher cumulants will approach equilibrium only after 
M and reach their respective equilibrium values. 


region. We note further that as inhomogeneities may be 
important for bubble nucleation in first order transitions, 
our studies will be restricted to the cross-over side of the 
critical point. 

In this work, we shall concentrate on the evolution of 
the fi rst fo ur cumulants. It is straightforward to extend 
Eqs. (2.20), if so desired, to include the evolution of even 
higher cumulants such as and kq. 


D. A brief summary 


We derived in this section a set of coupled first order 
differential equations, Eqs. (2.20), describing the evolu¬ 


tion of the first four cumulants, ft n ,n = 1,2, 3,4 of the 
zero mode a of the critical field in the vicinity of the 
critical point. A key feature of these equations is that 
the evolution of higher moments are only coupled to the 
evolution of lower moments, and not vice versa. One 
therefore obtains a closed form system of equations that 
can be solved numerically. 

In general, the evolution of the first four cumulants 
will be coupled to the evolution of higher cumulants, 
K n including n > 4 as well. However in the e <C 1 
limit, where are results are strictly applicable, we demon¬ 
strated analytically that the coupling of lower cumulants 
to higher cumulants is suppressed. Therefore the sys¬ 
tem of equations we derived is applicable for describing 
the temporal evolution of moments in the scaling regime 
Emicr < Ceq < E, where e C 1 is satisfied. For small 
values of e, we have checked explicitly that the difference 
between computing the non-equilibrium cumulants from 
numerical solutions of the Fokker-Planck master eq uation 
and from solutions of the evolution equations Eqs. (2.20) 
is suppressed by e (see Appendix. [B] for further discussion 
of this point). 

We could of course in pr inciple have solved the mas¬ 
ter equation in Eq. (2.10) directly. However, this is 


not advisable for both practical and conceptual reasons. 
Firstl y, fro m a practical perspective, numerically solving 
Eqs. (2.20) is much faster than solving the Fokker-Planck 


equation. Furthermore, the interplay between cumulants 
is far more transparent in the former approach. From 
a conceptual point of view, solving the Fokker-Planck 
equation Eq. (2.10) would not give a more faithful repre¬ 


sentation of how cumulants evolve in the crit ical r egime. 
This is because an important input into Eq. (2.10) is the 


equilibrium distribution function Po(cr), which is not easy 
to obtain. From universality, this distribution is also the 
equilibrium distribution of the 3-D Ising model. How¬ 
ever, while the cumulants in the 3-D Ising model are 
known, reconstructing Pq(ct) is non-trivial when e is not 
small. 

Thus conceptually there is no advantage in solving the 
Fokker-Planck equation at large e and no practical ad¬ 
vantage in solving it for small e. Hence for the current 
state of the art, the derived Eqs. (2.20) provide the most 


complete and consistent information, albeit limited, on 
the evolution of non-equilibrium moments in the critical 


III. OUT OF EQUILIBRIUM EVOLUTION OF 
CUMULANTS IN THE QCD CRITICAL REGIME 


In this section, we will apply Eqs. (2.20) to study the 
evolution of cumulants in the critical scaling regime of 
QCD. Our discussion is organized as follows. We will 
first discuss the problem in the context of the 3-D Ising 
model since it lies in the static universality class of the 
QCD critical point. This is important for determining 
the equilibrium values of the cumulants that provide the 
initial conditions for the evolution equations. In this con¬ 
text, we will also discuss the mapping of the Ising scal¬ 
ing regime to QCD as well as the dynamical universality 
class governing the transport properties of the medium 
near the critical point. Next, for the purposes of our 
study, we will construct a simple model of the medium 
that mimics the expanding fireball formed in heavy ion 
collisions. Finally, within the framework of this simple 
model, we shall present our results for the temporal evo¬ 
lution of cumulants along a representative trajectory in 
the QCD critical scaling regime. 


A. Fixing parameters from universality 


The inputs required to solve Eqs. (2.20) along a tra¬ 


jectory passing through the QCD scaling regime include 
(To, , A 3 , A 4 and the effective relaxation time r e ff. As 

we will also discuss shortly, we will need to know how 
V, T change with r for the case of an expanding medium. 

As noted previously, explicit expressions in the 3- 
D Ising model for the equilibrium ft® q (r, h) as a func¬ 
tion of the Ising variables are given in Appendix [A] 
With these expressions, we can determine the param¬ 
eters (Jo(r, h ), m a (r, h ), A 3 (r, h ), A 4 (r, h) from Eq. (2.9). 


The dependence of r e ff on £ eq is universal and can be 
expressed in terms of the dynamical critical exponent z 


(3.1) 

Here r re 1 can be interpreted as the relaxation time at 
the outside edge of the critical region, defined by a min¬ 
imal correlation length, £ = £ m i n . As we discussed ear¬ 
lier, the work of Ref. (26] demonstrated that transport 
properties in the vicinity the QCD critical point are gov- 


Teff — "Tel 
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FIG. 1. (Color online) A sketch of the cross-over side of the 
scaling regime in the r, h plane. We define the scaling regime 
with the criterion £ min < £ eq < £max. To be specific, we take 
(max/£min = 3. The solid curve delineates the boundary of the 
critical regime. The dotted curves show where the equilibrium 
value of the kurtosis K changes sign. The trajectory A (see 
text) is shown with the black dashed vertical line. 


erned by the diffusion of the baryon density Because 
the chiral order parameter is not conserved and mixes 
with the baryon density, it does not influence the dy¬ 
namical universality class. Consequently, the dynamical 
universality class is that of the liquid-gas phase transi¬ 
tion, model H; in the classification schem e of Ref. m, 
this gives z = 3. It is clear from Eq. CED that the relax¬ 
ation of the critical mode to equilibrium is greatly slowed 
down as a spacetime trajectory in the system approaches 
the critical point. 

Our discussion thus far only relied on the static and dy¬ 
namical universal properties of QCD critical point. How¬ 
ever, as we are interested in the evolution of cumulants in 
the QCD critical regime, we also need to specify the map[^] 
between the Ising variables r, h to the QCD variables 
T, /i. In the Ising model, the critical point is at r = h = 0. 
There is a first order transition for h = 0, r < 0 and a 
cross-over for h = 0, r > 0. Therefore the r-axis is the 
direction tangential to the line of first order transition 
ending at T c . In contrast, in general, the h-axis will de¬ 
form after the map to the T, /i plane. How this occurs is 
not known; for simplicity, we will follow the prescription 
of previous studies (21 HU and assume that h is perpen¬ 
dicular to the r-axis. Specifically, we will assume that h 
is parallel to the T axis in the QCD phase diagram. We 


2 It is also argued in Ref. that the conservation (or not) of the 
isospin density does not influence considerations of universality 
because the isospin susceptibility, unlike the baryon susceptibil¬ 
ity, is finite at the critical point. 

3 This map is non-universal which is a significant source of system¬ 
atic uncertainty in quantitative dynamical studies in the vicinity 
the QCD critical point. 


then have following linear mapping relations, 

T - T c h_ _ r_ 

AT Ah’ A /I Ar ’ 1 ' ’ 

where AT, A /jl denote the width of the critical regime in 
the QCD phase diagram. The corresponding width of the 
critical regime in the Ising variables Ar, Ah is defined to 
be 


£(r = Ar, h = 0 ) = £(r = 0 , h = Ah) = £ n 


(3-3) 


The inner and outer boundary of the critical regime in the 
r, h plane are illust rate d in Fig. [T| w here the curves are 
obtained from Eq. (A5) and Eq. (A2), for £ m ax/£min = 3. 
The equilibrium kurtosis is even with respect to h. Due 
to the potential importance of observables for the sign 
of k urtos is, we als o pl ot in Fig. [l] the boundary (from 
Eq. (A5) and Eq. ( |A2| )) where the equilibrium kurtosis 
flips sign. It is worth pointing out that, in contrast, the 
equilibrium skew ness i s an odd func tion of h. One may 
check from Eqs. (|A2|), (|A3|) and Eq. (|A5|), that the skew¬ 


ness is positive for ne gativ e h. Given our choice of the 
direction of h -axis in (3.2), the equilibrium skewness is 


negative above the cross-over and is positive below it. 


B. Simple model of space-time evolution in the 
vicinity of the critical point 


For a given center of mass energy (yT), the space-time 
trajectory in the QCD phase diagram, can be determined 
by hydro dynamical simulations if the expansion rate is 
smaller than characteristic scattering rates. For the il¬ 
lustrative purposes of the study in this paper, a simple 
dynamical model of space-time evolution is sufficient. We 
will assume that for fixed yT, the /i of the fireball is con¬ 
stant during the evolution. This trajectory corresponds 
to the vertical dotted line in the critical regime shown 
in Fig. [Tl It implicitly assumes of course the mapping 
relation Eq. (O). 


We parametrize the evolution of the volume as 


V(r) 

Vt 


Tl 


T(t) 

T t 


Tl 


(3.4) 


where Vi and Tj are respectively the volume and tem¬ 
perature of the system at tj, the time at which the tra¬ 
jectory first hits the boundary of critical regime. Here 
ny controls the rate of expansion; ny = 3 corresponds 
to a 3 dimensional Hubble-like expansion and ny = 1, 
a 1 dimensional Bjorken-like expansion. To obtain T(r) 
in Eq. (3.4), we assumed that the total entropy is ap¬ 
proximately conserved during the evolution, and hence 
the entropy density goes as s(r) ~ r~ nv . We also used 
the relation dlogT/dlog s = c 2 3 , where c 2 is the speed 
of sound. Even at the highest heavy ion collision ener¬ 
gies, the system is more Hubble-like than Bjorken-like 
when trajectories approach the critical point. We will 






















therefore take ny = 3 for our study. Further, guided 
by lattice measurements [39h42l which indicate that c 2 s is 
around 0.15 near the QCD cross-over line, we will pick 
this value for the study in this paper. 

Finally, we need to sp ecify the initial conditions for 
the solution of Eqs. (2.20). We will assume initially that 
M, ft n , n = 2, 3,4 are equal to their equilibrium values at 
r = rj. Choosing the initial volume is tricky, because, for 
the reasons outlined previously, we wish to ensure that 
e< 1 at all times. We do this by requiring that the max¬ 
imal value of e is imposed to be e c e 


max 

v c 


= 0.1. Here 

V c is the volume when the trajectory followed by the sys¬ 
tem crosses T c in the cross-over region. As long as the 
change of the equilibrium correlation length is faster than 
the expansion of the volume, e c gives the upper bound 
of e during the evolution of the medium^] For instance, 
for Cmax = 3 fm, e c = 0.1 corresponds to V c ~ (14 fm ) 3 . 
We ca n the n use the time evolution of the temperature 
in Eq. (3.4) to determine tj, and subsequently the corre¬ 
sponding equation for V (r) to determine Vi in terms of 
ti- 

With these mode l assumptions, we are now ready 
to solve Eqs. (2 .20| ) for each given trajectory passing 
through the critical regime. Our results will only de¬ 
pend on one dimensionless parameter, r re \/rj , where r re \ 
is the relaxation time at the boundary of critical regime. 
As there are no extant first principles calculations of r re i 
in the QCD critical regime, w e sha ll take its value as a 
free parameter and solve Eq. (2.20) for different choices 
of r re i/r/. We can benchmark this value by noting that 
if i) r re i is 1 fm, a characteristic strong interaction scale, 
and ii) r/, the initial time at which the system enters the 
critical regime, is 10 fm, a reasonable estimate would be 

Trel/T/ ~ 0.1. 


C. Results for cumulant evolution along a 
representative trajectory 


For later convenience, we will express the evolution of 
the non-Gaussian cumulants in terms of the dimension¬ 
less quantities skewness (S) and kurtosis (K): 


S = 


fVl/2 

V 4 k 3 

3/2 


K = 


V 4 K 4 


(3.5) 


where V 4 is the rescaled V 4 , 


(3 - 6) 

With these definitions, we have deliberately taken the 
four-volume dependence out, as one should in explo¬ 
rations of critical behavior. 


4 For e< 1, our results are independent of the choice of e. 


We will first study the evolution of the non-equilibrium 
cumulants along the representative trajectory A in Fig. [l] 
Along this particular trajectory, £ eq will approach £ max 
when the trajectory approaches the cross-over line. Mem¬ 
ory effects are therefore most prominent along trajectory 
A. For each fixed /i, trajectories can be parametrized by 
h, or equivalently T, via the mapping previously out¬ 
lined. We plot in Fig. [2| as a function of varying T, 
the non-equilibrium ratios M/Ma, £/£min, S/Sa, K/Ka- 
Here Ma , £ m in , Sa , Ka are the equilibrium values speci¬ 
fied at the end point of the trajectory. 

One immediately observes that non-equilibrium effects 
are important for all cumulants. The difference be¬ 
tween the non-equilibrium cumulants and equilibrium 
cumulants is apparent unless t tg \/ti is much smaller 
than the noted benchmark value. (See for instance the 
red curves in Fig. [2] that correspond to solutions with 
r re i/ t i = 0.005.) The difference is most visible near the 
cross-over line (T = T c ) where £ eq reaches its maximum 
value along the trajectory. This deviation is a direct 
manifestation of the effects of critical slowing down. 


We shall now examine in detail the evolution for each 
individual non-equilibrium cumulant. The behavior of 
M(T) is relatively simple. As Fig. [ 2 ] (a) shows, M 
tends to follow the equilibrium values, though with crit¬ 
ical slowing down, the change in sign occurs later for 
larger values of r re 1 . £ eq is an even function of the Ising 
variable h-due to the mapping relation Eq. (3.2)—it is 
symmetric with respect to T c . Our results for the non¬ 
equilibrium correlation length are qualitatively similar 
to p reviou s studies based on solving the rate equation 
Eq. (2.23) [24[ [27] . The role of memory effects on £ is 
two-fold. On the one hand, when the equilibrium correla¬ 
tion length £ eq is large, the effects of critical slowing down 
delay the growth of the effective non-equilibrium correla¬ 
tion length £. For example, as shown in Fig. [ 2 ] (b), when 
T is close to T c and £ eq approaches its maximum, the non¬ 
equilibrium £ for all r re 1 under consideration are smaller 
than the equilibrium value. On the other hand, memory 
effects of the critical regime are preserved more efficiently 
than if the system were in equilibrium throughout. One 
observes that when T is below T c , the non-equilibrium 
value of £ at that temperature is larger than the equilib¬ 
rium value. Similar observations were made previously 
in Ref. l24l. 


Turning now to the evolution of the non-Gaussian cu¬ 
mulants S and K, we first recall that in equilibrium /Cg q , 
or equivalently S eq , is an odd function of the Ising vari¬ 
able h. It will therefore will flip sign when crossing the 
cross-over line, as demonstrated by the dashed curve in 
Fig. [ 2 ] (c)). In contrast, K e 4 q or K eq is an even function 
of the Ising variable h. It is negative at the cross-over 
temperature and positive away from it as shown in the 
corresponding dashed curve in Fig. [ 2 ] (d). However we 
demonstrate in Fig. [2] that the non-equilibrium evolution 
of skewness and kurtosis do not necessarily follow the evo¬ 
lution of the corresponding equilibrium cumulants, and 
can be radically different in both their magnitude and 
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FIG. 2. (Color online) (From top to bottom) The evolu¬ 
tion of non-equilibrium mean M/Ma{ a), effective correlation 
length £/£min (b), skewness S/Sa (c) and kurtosis K/Ka (d) 
as a function of (T c — T)/AT along trajectory A (c.f. Fig.[l]). 
Results for T re i/rj = 0.005,0.02,0.05,0.2 are shown in solid 
red, dotted blue, single dot dashed green, double dot dashed 
orange curves respectively. The dashed curves plot the cor¬ 
responding equilibrium values. All results are normalized by 
the corresponding equilibrium value at the end point trajec¬ 
tory A (c.f. Fig. [I]). The dashed vertical lines illustrate the 
T at the point where the a freeze-out curve intersects with 
trajectory A (from left to right, freeze-out curves o f type I, 
II, III respectively. For a discussion, see text in Sec. IV B ). 


sign. 

These differences occur because, as previously noted, 
the evolution of higher cumulants is coupled to the lower 
ones. Therefore how K (or S) evolve will not only de¬ 
pend on its deviation from the equilibrium value, but also 
on the non-equilibrium values of other cumulants. As 
we shall discuss shortly, these deviations off-equilibrium 
have significant phenomenological implications for the 
search for a critical point in a beam energy scan. Specif¬ 
ically, in Fig. [2| the dashed vertical lines correspond to 
freeze-out trajectories I, II, and III (left to right), which 
provide snapshots of the non-equilibrium cumulants that 
may be measured in experiments. We shall return to a 
more detailed discussion of these in Sec. WE 


IV. TOWARDS MODELING THE RHIC BEAM 
ENERGY SCAN 

The results we presented for the non-Gaussian cumu¬ 
lants off-equilibrium potentially strongly impact the in¬ 
terpretation of the results of ongoing and future criti¬ 
cal point searches with the beam energy scan (BES) at 
RHIC. To further explore these, we will solve the evolu¬ 
tion equation for fixed-// trajectories broadly spanning 
the critical regime. In our simple model, this would 
be the equivalent of varying y/s. We will then be able 
to compute the non-equilibrium cumulants for a given 
r re i/ t i for every point in the critical regime. 


A. Memory effects and the sign of 
non-equilibrium skewness and kurtosis 

In Figs. [3] and [4] respectively, we present contour plots 
for the equilibrium and non-equilibrium skewness S and 
kurtosis K. Due to memory effects, the non-equilibrium 
contours in the T — fi plane deform from the correspond¬ 
ing equilibrium contours; the deformation is enhanced for 
larger relaxation times r re i /tj . 

We now focus on the sign of the skewness and kurtosis, 
their most prominent feature. We illustrate it by plotting 
regime S > 0 (or K > 0) in red and S < 0 (or K < 0) 
in blue. In equilibrium, the boundary that separates the 
regime where S > 0 and S < 0 is precisely the cross-over 
line at T = T c . In Fig. [3j we fix the sign of the equilib¬ 
rium skewness in such a way that S eq is positive below 
the cross-over line; this is suggested by the arguments 
presented in Refs. mm and the observation that the 
skewness is positive in the hadron resonance gas model. 
Fig. [3] demonstrates that for non-equilibrium skewness S 
the boundary where the skewness changes sign deforms 
and with increasing r re \/rj becomes negative in a larger 
portion of the area below the cross-over line. This is to 
be expected as the non-equilibrium cumulants carry more 
information at early times when the equilibrium skewness 
is negative; a larger relaxation time r re \ would give more 
weight to early time contributions. Similarly, as Fig. [4] 
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FIG. 3. (Color online) Contour plot of equilibrium (top) and 
non-equilibrium skewness (S) with r/rj = 0.05 (middle) and 
t/tj = 0.02 (bottom). The S > 0 region is shown in red and 
S < 0 region is shown in blue. 


FIG. 4. (Color online) Contour plot of equilibrium (top) and 
non-equilibrium kurtosis (K) with r/rj = 0.05 (middle) and 
r /tj = 0.02 (bottom). The X > 0 region is shown in red and 
K < 0 region is shown in blue. 


shows, the boundary separating the regime K > 0 and 
K < 0 also deforms. 

With Figs. [3] and [4] in mind, we may ask how the sign 
of the non-equilibrium skewness and kurtosis behaves as 
a function of /jl (or y/s) on the freeze-out curve of a heavy 
ion collision in the T-fi plane. As Fig. [3] illustrates, the 
skewness on the freeze-out curves can be either negative 


or positive depending on the magnitude of r re i/r/ and 
the relative position of the freeze-out curves. However, 
despite memory effects, the sign of non-equilibrium kur¬ 
tosis will still switch from negative to positive when /i 
approaches /i c along the direction of cross-over. This 
change in sign of the /1 dependence of kurtosis is insensi¬ 
tive both to the magnitude of r re 1 as well as the location 
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FIG. 5. A sketch of the cross-over side of the scaling regime in 
the T, fi plane. The dotted curves show where the equilibrium 
value of the kurtosis K changes sign. The trajectory A (see 
text) is shown with the black dashed vertical line. Different 
types of freeze-out curves (F.C.) are shown in red (upper ), 
blue (middle) and green (lower) dashed curves, corresponding 
to type I, II, III respectively (see text). 


of the freeze-out curves. 


B. Dependence of non-equilibrium cumulants on 
fs and freeze-out curves 


We shall now illustrate a few possible experimental 
outcomes for the behavior of cumulants as a function 
of fi (or fs) if the cross-over side of critical regime of 
QCD phase diagram is scanned. We emphasized pre¬ 
viously that the cumulants of the critical field a itself 
are not directly observable. However the critical field a 
is coupled to the net baryon number density and there¬ 
fore critical fluctuations contribute to the moments of 
net baryon number fluctuations that are measured in ex¬ 
periments. Indeed, such contributions are proportional 
to the corresponding moments of the cr field itself mm- 
We may therefore expect that the //-dependence of cumu¬ 
lants of the critical fields, as determined in our model, will 
capture the qualitative behavior of cumulants of particle 
fluctuations. 

Consider Fig. [5j where we have now re-plotted Fig. [T| 
perfo rmin g the map from r, h variables to T and // using 
Eq. (3.2). In this figure, we have superposed the tra¬ 
jectories corresponding to three freeze-out curves. Each 
different choice of a freeze-out curve corresponds to tak¬ 
ing a different snapshot of the evolution of cumulants as 
represented by the dashed vertical curves in Fig.[2|. Our 
results will therefore depend on the relative overlap be¬ 
tween the QCD critical regime and the particular freeze- 
out curve in the QCD phase diagram. To describe the 
freeze-out curves, we take an empirical parametrization 
of the heavy-ion collision data from Ref. [43] , 


7 >( m ) = a - b/i 2 B - c/j.% , 


(4.1) 


with a = 0.166 GeV, b = 0.139 GeV -1 , c = 0.053 GeV -3 . 


Given the mapping (3.2), the overlap between the crit¬ 
ical regime and the freeze-out curves depends on // C ,T C 
as well as AT, A//. Currently the location of the QCD 
critical point and the width of QCD critical regime are 
not known. Model calculations [44] and lattice QCD cal¬ 
culations [45! suggest that A/i ~ 0.1 GeV. We therefore 
set ji c = 0.25 GeV, A/i = 0.1 GeV, A T/T c = 1/8. Conse¬ 
quently, the overlap between a freeze-out curve and QCD 
critical regime will depend on T c . In practice, we shall 
take three different values, T c = 0.165,0.18,0.194 GeV, 
to represent three freeze-out trajectories overlapping with 
the critical regime: I) the freeze-out curve is near the 
cross-over line, II) below the cross-over line and III) near 
the edge of the critical regime. They are plotted in Fig. [5] 
and will be labeled I, II, III respective!^] 

For each constant /i, we will denote the values of the 
non-equilibrium cumulants at the point where the tra¬ 
jectory intersects the freeze-out curves by the subscript 
“F”. We will then study the dependence of the cumulants 
on freeze-out curves as a function of fi. Our results for 
£f(/x), SV(aO, Kp{f) are shown in Figs. [6] [7] [8] respec¬ 
tively. We observe generally that the // dependence of 
the non-equilibrium cumulants can be different from the 
equilibrium cumulants (which are represented by dashed 
lines). As we anticipated previously, the values of the 
non-equilibrium £ on all the freeze-out curves are consid¬ 
erably amplified for a wide range of r re i/r/. This implies 
that even if the freeze-out curve is located at the edge of 
the critical regime (as for the freeze-out curve III shown 
in Fig. |6] (c)), memory effects ensure that the signature 
of critical fluctuations is not necessarily suppressed. Re¬ 
garding skewness, we noted arguments advanced that the 
equilibrium skewness is positive below the cross-over line. 
Our results in Fig. [7] clear demonstrate that memory ef¬ 
fects can modify the sign of skewness off-equilibrium to 
be opposite to that of equilibrium skewness. Similar de¬ 
viations from equilibrium expectations are observed for 
kurtosis, as shown in Fig. [8] 

From Figs. [7] and [8j we observe that the non-Gaussian 
cumulants Sf{l0 and Kpf) are sensitive to the relative 
positions of freeze-out curves (or T c ) and r re We can 
go one step further, and given a model for the fs de¬ 
pendence of //, examine how these fluctuations vary as 
a function of fs. Ref. [43], from fits to data, developed 
the parametrization, 


li(Vs) = 


do 

d\fs + 1 


(4.2) 


5 A reader might question whether the values of T c considered 
are reasonable given our present knowledge of the QCD phase 
diagram. However, for our illustrative purposes, the choice of 
T c is just an easy way to represent the distance of the critical 
point from the freeze-out trajectory. Alternative ways of rep¬ 
resenting this relative separation are certainly feasible; for in¬ 
stance, we could have chosen smaller values of T c and modified 
the parametrization of the freeze-out curves. 
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FIG. 6. (Color online) The non-equilibrium values of the ef¬ 
fective correlation length on the freeze-curves as a function 
of fi. The results corresponding to freeze-out curves of type I, 
II, III (as shown in Fig. [5]) are (a), (b), (c) respectively. Non- 
equilibrium cumulants with T re i/r/ = 0.005,0.02,0.05,0.2 are 
represented by solid red, dotted blue, single dot dashed green, 
double dot dashed orange curves respectively. The dashed 
curves plot the corresponding equilibrium values-where not 
visible, they fully overlap with the solid red curve. All results 
are normalized by the corresponding equilibrium value at the 
end point of trajectory A as shown in Fig. [5] 
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FIG. 7. (Color online) The non-equilibrium values of the 
non-Gaussian cumulant Sf on the freeze-curves as a function 
of fi. The results corresponding to freeze-out curves of type I, 
II, III (as shown in Fig. [5]) are (a), (b), (c) respectively. Non- 
equilibrium cumulant results for T re \/ri — 0.005, 0.02, 0.05, 0.2 
are displayed with solid red, dotted blue, single dot dashed 
green, double dot dashed orange curves respectively. The 
dashed curves plot the corresponding equilibrium values. All 
results are normalized by the corresponding equilibrium value 
at the end point trajectory A, as shown in Fig. [5}. 


where do = 1.308 GeV and di = 0.273 GeV -1 . Using 
this relation, we obtain the results shown in Fig. [9] for 
the skewness and kurtosis respectively. In particular, we 
demonstrate that very similar curves, as a function of y ds 
can be obtained by different combinations of freeze-out 
curves and relaxation times. These results suggest that 
great care should be exercised in interpreting trends as 
a function of y/s of measured non-Gaussian cumulants. 
While these may signify the onset of critical dynamics, 
more needs to be done to model freeze-out conditions and 


constrain r re i before definitive statements can be made 
regarding the discovery of a critical point. 


SUMMARY AND CONCLUSION 


We derived in this paper a set of equations Eqs. (2.20) 


that describe the evolution of non-equilibrium cumu¬ 
lants of the critical field a in the QCD critical regime. 
In particular, we obtained novel expressions for the 
off-equilibrium evolution of non-Gaussian cumulants. 
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FIG. 8. (Color online) The non-equilibrium values of the 
non-Gaussian cumulant Ff on the freeze-curves as a function 
of /i. The results corresponding to freeze-out curves of type I, 
II, III (as shown in Fig. [5]) are (a), (b), (c) respectively. Non- 
equilibrium cumulant results for r re i = 0.005, 0.02, 0.05, 0.2t/ 
are displayed with solid red, dotted blue, single dot dashed 
green, double dot dashed orange curves respectively. The 
dashes curves plot the corresponding equilibrium values. All 
results are normalized by the corresponding equilibrium value 
at the end point trajectory A, as shown in Fig. [5|. 


With equilibrium initial conditions given by the three- 
dimensional Ising model (which belongs to the static 
universality class of the QCD critical point), we applied 
these equations to study the real time dynamics of cu¬ 
mulant s along a trajectory passing through the cross-over 
side of the critical regime. Since this is a first exploratory 
study, our studies were performed in the framework of 



18 20 22 24 26 


Vs (GeV) 


FIG. 9. (Color online) The yfs -dependence (assuming 
Eq. (4.2)) of non-equilibrium non-Gaussian cumulants on the 
freeze-curves with different choices of the relative position of 
freeze-out curves and r re i- (Top) skewness & vs y/s; (Bot¬ 
tom) kurtosis Kf v s y/s. The results obtained from the com¬ 
bination, i.e., (F.C. I, T//r re i = 0.02), (F.C. II, rj/r re i — 0.05), 
(F.C. Ill, Ti/r re i = 0.2), are displayed in solid red, dotted 
blue, single dot dashed green curves respectively. 


a very simple model of expansion dynamics in heavy- 
ion collisions. With further data motivated parametriza- 
tions of chemical freeze-out in the T — (i plane, we were 
able to obtain snapshots of the non-equilibrium cumu¬ 
lants along freeze-out curves that overlapped to differing 
extents with the QCD critical region. 

Our chief conclusion is that memory effects are respon¬ 
sible for substantial differences in the temporal evolution 
of non-equilibrium and equilibrium cumulants. Specifi¬ 
cally we found that the sign of non-equilibrium skewness 
in the vicinity of the critical point can be opposite to 
the value expected in equilibrium. Similar conclusions 
are suggested by our results for non-equilibrium kurto¬ 
sis. Our results are sensitive to the relaxation time r re i 
of the critical field a fluctuations and the location of the 
freeze-out curves in the QCD phase diagram relative to 
the extent of the QCD critical region. We find that mem¬ 
ory effects are important to such an extent that even tra¬ 
jectories that traverse the edge of the critical regime are 
sensitive to its dynamics. 

Our treatment of critical fluctuations in the QCD crit¬ 
ical regime was strongly motivated by the description of 
similar fluctuations in the three-dimensional Ising model. 
In the latter case, cumulants of critical fluctuations are 
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expressed in terms of the reduced temperature r and the 
rescaled magnetic field ft. The map of the description 
of critical fluctuations in terms of r and ft to T and (i 
is non-universal, and is a significant source of system¬ 
atic uncertainty in treatments of critical dynamics in the 
QCD critical regime. This uncertainty, coupled with our 
ignorance of r re i, provide fundamental obstacles to quan¬ 
titative studies of real time critical dynamics in QCD. 

Indeed, because of the importance of non-equilibrium 
effects, lattice studies of equilibrium cumulants, while of 
fundamental importance, may not be sufficient. These 
must be accompanied by progress in non-equilibrium 
studies of the QCD critical regime. One promising ap¬ 
proach is the use of classical statistical real time simula¬ 
tions mu mi that have also previously been applied to 
studying the non-equilibrium dynamics of the very ear¬ 
liest stages of high energy heavy-ion collisions 03 ng. 
Detailed dynamical models of the space-time evolution 
of heavy-ion collisions as a function of beam energy are 
also very important. In particular, models that build in 
the transport of conserved charges and reproduce bulk 
features of these collisions such as particle spectra can 
place strong constraints on the parameter space for the 
non-equilibrium evolution of cumulants. 

In this work, we have concentrated on critical dynam¬ 
ics on the cross-over side of the critical regime. From the 
perspective of a critical point search, this approach is ap¬ 
propriate because it is easier in both experiments and in 
lattice gauge theory computations to extend explorations 
of the QCD phase diagram starting from the regime of 
high temperatures and low baryon chemical potentials. 
However, if a critical point is localized, it would be of 
great interest to understand non-equilibrium dynamics 
on the first-order side of phase diagram. In this regard, 
applying the framework discussed here from the cross¬ 
over critical regime to the first-order critical regime of 
the QCD phase diagram is a useful extension to be pur¬ 
sued in future studies. 
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rium cumulants can be computed by taking derivatives 
of M eq (r, ft) with respect to ft at fixed r, 


eq = 1 fd n M^(r,h)\ 

' n+1 (V 4 H 0 ) n V dh n ) r ' 


n = 1,2,3,... 


(A!) 

Here Hq is a dimensionful parameter (of mass dimen¬ 
sion 3) which relates reduced magnetic field ft to the un¬ 
reduced magnetic field. 

To parametrize M eq (r, ft), we use the linear parametric 
model EH3 ED]. In this parametrization, one introduces 
two new variables 77, 6 which are related to (dimension¬ 
less) Ising variable r, ft as 


r(R, 9) = ,R(1 — 6 2 ), h(R, 6 ) = Ah R 136 h{9 ), (A2) 


Following Ref. m , we will use 


h{0) = 30 



{ («-3) ) S 


(A3) 


Here /?, S are standard critical exponents and we will use 
the values obtained from mean field theory, (3 = 1/3, 5 = 
5. In these 77, 6 variables, 0 = 0 corresponds to the 
crossover line and \ 6 \ = a/3/2 corresponds to the co¬ 
existence (first order transition) line. The equilibrium 
“magnetization” Mg q (r, ft) (or ctq) is given by 


M eq (77, 6 ) = MqR^O , (A4) 


where M 0 sets the scale of “magnetization”. The 
parametrization introduced describes the equation of 
state with a precision sufficient for our purpose. 


We now compute K® q using Eq. (|A1| and Eq. (A4). 
Explicitly, we have 


/4W) 

/4W) 

/COM) 


Mo 1 

R 4 M 0 R 4 /3(3 + 26> 2 ) ’ 

-M 0 40(9 + 6 2 ) 

(V 4 H 0) 2 R 3 (3 — 6 2 )(3 + 2# 2 ) 3 ’ 
— 12M 0 

(W 


(A5) 

(A6) 


X 


(81 - 783<9 2 + 105<9 4 - 5<9 6 + 2<9 8 ) 

Ri 4 /3(3_ #2)3(3 + 2£ 2 ) 5 


.(A7) 


Fina lly, we convert /s® q (i7,0) into ft® q (r, h) using 
Eq. (A2). We note that M/Ma, C/Cmin, S/Sa, K/Ka as 
presented in this paper does not depend on the choice of 
dimensionful normalization Mq , Hq . 


Appendix A: Parametric representation of 
equilibrium cumulants in the Ising critical regime 


Appendix B: Detailed derivation of Eqs. (2.20) 


In this section, we explain the parameterization of the 
equilibrium cumulants M eq (r, ft), ft® q (r, ft), n = 2, 3,4,... 
in the critical regime in terms of the Ising variables r and 
ft used in this paper. For this purpose, we only need to 
know the equilibrium magnetization M eq (r, ft) as equilib- 


We present here a detailed derivation of Eqs. (2.20). 
It is convenient to introduce the generating function of 
cumulants, 


G( A; r) = log [Z{ A; r)] . Z( A; r) s= (e A ^>. (Bl) 
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where the average (...) and 5 a have been defined in 


given by 


Eq. (2.2) and Eq. (|2.3|) respectively. The cumulants are 


<9 t k 2 (t) = -2 rj(b 2 ) 


- 1 




d n G{ A) 


d\ n 





e 2 

+ V 

|_Ve 6 2 / 

? 

n = 2,3,4,... . 

(B2) 

A 


A=0 



+ 6 1 

f «4 \ 




+ _ 6 1 

Ke 2 b 4 ) 


We also have 5aG(A;t)|a=o = 0 as (5a) = 0. Equiva¬ 
lently, one could read ft n by Taylor expanding G(A; r) 
around A = 0 , 


G(X;t) = ( B3 ) 

n=2 

We now take the derivative of G'(A; r) with respect to 


d T G(\:r) = 


8 t Z(X;t) 

Z(\;t) 


= E • 


(B4) 


n—2 


Therefore the evolution equations for ft n can be deter¬ 
mined by evaluating 3 r G(A;r). We consider 

d T Z(X; t ) = d T [e~ XM (e Xa )] 

= -AZ(A; t)(8 t M) + e~ XM d T {e Xa ) . (B5) 


A 2 


From Eq. (2.11), we have 
e~ XM d T (e Xa ) =-T_ 

m a T eff 

(B 6 ) 

We can now substitute Eq. (B 6 ) into Eq. (B5) and then 
plug the results into Eq. (B4) to obtain, 


X(e x ^%(a)) - -Z(X;r) 


8 t G(X; t) = -Xd T M - — J F,(M) 

'"'Tctr 


e(5cre XScr \ 


bZ 

f e 3 (5<j 3 e XScT ) 


F 2 (M) 


2 (Sc 


b 2 Z 


— Xb 


Fs(M) 

(B7) 


Now the evolutions equations for M and cumulants 
K n , n = 2,3,4,... can be determined by Taylor expand¬ 
ing both sides of the above equation in powers of A 
and comparing coefficients in front of A". At order 
A", n = 1 , 2 ,3,4, we have, respectively, 


8 t M(t) 



? (5) F ‘ (M) ■ 


Fi(M) + e -^)F s (M) 

(B 8 a) 


8 t k 3 (t) : 


+ 


-3r eff 1 (eb 3 )| 




^4(t) = -4^ <j [(^j)F 2 (M) + 3 (^)i^)Fs(M) 


f 1^2 \ / ^3 


+ (>. 


e 2 

+ — 


(5^3 (M) 


2 LW 


+ (»<£>*+ <><&> I* 


(B 8 d) 


In deriving these equations, we also used the relations 

(5a 2 e XSa ) 


( Sae XScr ) 


= 3 A G, 


— \p\G + (8\G) 2 


= [8lG + 3(d 2 G)(d x G) + (8 x G) s 


(B9) 


(5a 3 e XSa ) 

~Z 

and Eq. 

Keeping contrib ution s to leading order in e in Eq. (|B 8 j) , 
we arrive at Eqs. (2.20). 

We note that Eqs. (B 8 ) would still be closed if we 


further include e 2 terms in the evolution equations for 


M, ft 2,^3 (but neglect e 2 terms in Eq. (B 8 d)). The re¬ 


sulting equations, which we shall refer as the next to 
leading order (NLO) evolution equations, take full sub¬ 
leading contributions in e 2 to the evolutions of M, ft 2 , ft 3 
and part of sub-leading contributions in e 2 to the evolu¬ 
tion of ft 4 into account, neglecting the term proportional 
to fts. We have checked numerically that the differences 
in computing the non-equilibrium cumulants from solu¬ 
tions of the Fokker-Planck master equation Eq. ( | 2 T 0 | , 
relative to i) sol utions of the leading order (LO) evolution 
equations Eqs. (2.20) and ii) from NLO evolution equa¬ 
tions become smaller and smaller with decreasing e 2 . For 
fixed 6 , the difference is relatively larger for non-Gaussian 
cumulants. In particular, since the sign of non-Gaussian 
cumulants ft 3 ,ft 4 is in-definite, we observed numer ically 
that occasionally, ft 3 ,ft 4 as determined from Eqs. (2.20) 


would oscillate around zero for a short period of r. In¬ 
deed, such behavior will disappear if one solves NLO evo¬ 
lution equations, which stabilizes the solutions against 
these oscillations. For this reason, the results presented 
in this paper are determined in practice by solving NLO 
evolution equations. 
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